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A parallel processing architecture and method, based on the 
3D Algebra Reconstruction Technique (ART), for iteratively recon- 
structing from cone beam projection data a 3D voxel data set (V) re- 
presenting an image of an object. The cone beam projection data is 
acquired by employing a cone beam x-ray source and a two-dimen- 
sional array detector to scan the object along a source scanning trajec- 
tory to obtain, at each of a plurality of discrete source positions (0i) 
on the source scanning trajectory, a 2D measured cone beam pixel 
data set The 3D voxel data set (V) is subdivided into or organ- 
ized as a plurality (m) of independent voxel subcubes (VO) through 
(Vm-i) each containing a plurality of voxels. As a result of the sub- 
division of the 3D voxel data set (VHnto voxel subcubes, the 2D 
measured cone beam pixel data set (^) (measured projection data 
array) is correspondingly subdivided for each source position (0,) 
into projection data subsets, with overlapping regions between one 
or more adjacent projection data subsets corresponding to adja- 
cent voxel subcubes. Each voxel subcube and its corresponding 
projection data strip or subset is processed in parallel with other 
pairs of voxel subcubes and corresponding projection data strips 
or subsets, without interference. A bi-level parallel-processor archi- 
tecture is employed for iterative reproject and merge, and split and 
backproject operations. 
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PARAIXEL PROCESSING METHOD AND APPARATUS BASED ON TEffi ALGEBRA REOONSTRUCTION 
TBCHNiaOB FOR RECXXISIROCIING A ![HREE-DINENS10NAL OQHPUIERIZED TOMOGE^APHY 



Background of thfi Invention 

The present invention relates generally to three- 
dimensional (3D) computerized tomography (CT) and, more par- 
ticularly, to methods and apparatus employing parallel pro- 
cessing for reconstructing a 3D image of an object from cone 
IC beam projection data. 

In conventional computerized tomography for both 
medical and industrial applications, an x-ray fan beam and a 
linear array detector are employed. Two-dimensional (20) 
imaging is achieved. While the data set is complete and 

15 image quality is correspondingly high, only a single slice of 
an object is imaged at a time. When a 3D image is required, 
a "stack of slices" approach is employed. Acquiring a 3D 
data set a 2D slice at a time is inherently tedious and time- 
consuming. Moreover, in medical applications, motion arti- 

20 facts occur because adjacent slices are not imaged simultane- 
ously. Also, dose utilization is less than optimal, because 
the distance between slices is typically less than the x-ray 
collimator aperture, resulting in double exposure to many 
parts of the body. 

25 A more recent approach, based on what is called 

cone beam geometry, employs a two-dimensional array detector 
instead of a linear array detector, and a cone beam x-ray 
source instead of a fan beam x-ray source, for much faster 
data acquisition. To acquire cone beam projection data, an 

30 object is scanned, preferably over a 360* angular range, 

either by moving the x-ray source in a circular scanning tra- 
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jectory, for example, around the object, while keeping the 2D 
array detector fixed with reference to the source, or by 
rotating the object while the source and detector remain sta- 
tionary. In either case, it is relative movement between the 
5 source and object which effects scanning, compared to the 
conventional 2D "stack of slices" approach to achieve 3D 
imaging of both medical and industrial objects, with improved 
dose utilization. 

However, image reconstruction becomes complicated 
10 and requires massive computational power when a 3D image is 
reconstructed from cone beam projection data, in contrast to 
reconstruction of a 2D image from fan beam projection data. 

Literature discussing the cone beam geometry *for 3D 
imaging and generally relevant to the subject matter of the 

15 invention includes the following: Gerald N. Minerbo, "Convo- 
lutional Reconstruction from Cone-Beam Projection Data", IEEE 
Trans. Nucl. Sci., Vol. NS-26, No. 2, pp. 2682-2684 (April 
1979); Heang K. Tuy, "An Inversion Formula for Cone-Beam 
Reconstruction", SIAM J. Math., Vol. 43, No. 3, pp. 546-552 

20 {June 1983); L.A. Feldkamp, L.C. Davis, and J.W. Kress, 

"Practical Cone-Beam Algorithm", j. Opt. Soc. Am. A., Vol, 1, 
No. 6, pp. 612-619 (June 1984); Bruce D. Smith, "Image 
Reconstruction from Cone-Beam Projections: Necessary and 
Sufficient Conditions and Reconstruction Methods", IEEE 

25 Trans. Med. Imag., Vol. MI-44, pp. 14-25 (March 1985); and 

Hui Hu, Robert A. Kruger and Grant T. Gullberg, "Quantitative 
Cone-Beam Construction", SPIE Medical Imaging III: Image 
Processing, Vol. 1092, pp. 492-501 (1989). In general, this 
literature discloses various formulas for reconstruction of 

30 an image, including the use of a 3D Fourier transform or a 3D 
Radon transform. Questions of data completeness achieved 
with various source scanning trajectories are also con- 
sidered. 
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The present invention is an implementation of the 
Algebra Reconstruction Technique (ART) for reconstructing an 
image of an object from its projections. The ART is an iter- 
ative algorithm which uses raysums to correct each voxel 
5 (volume element) of the reconstructed image at each itera- 
tion. The ART was initially disclosed in Richard Gordon, 
Robert Bender and Gabor T. Hemnan/ "Algebraic Reconstruction 
Techniques (ART) for Three-dimensional Electron Microscopy 
and X-ray Photography", J. Theor. Biol., 29, pp, 471-481 

10 (1970) . Details regarding the mathematical foundations of 
the ART algorithm were described in Gabor T. Herman, Arnold 
Lent and Stuart Rowland, "ART: Mathematics and Applica- 
tions — A Report on the Mathematical Foundations and on the 
Applicability to Real Data of the Algebraic Reconstruction 

15 Techniques", J. Theor. Biol., 43, pp. 1-32 (1973). Extension 
of the additive ART to 3D cone beam geometry is described in 
M. Schlindwein, "Iterative Three-Dimensional Reconstruction 
from Twin-Cone Beam Projections", IEEE Trans, Nucl. Sci . , 
Vol. NS-25, No. 5, pp. 1135-1143 (October 1978). The ART was 

20 recently used for vascular reconstruction, as reported in A. 
Rougea, K.M. Hanson and D. Saint -Felix, "Comparison of 3-D 
Tomographic Algorithms for Vascular Reconstruction", SPIE, 
Vol. 914 Medical Imaging II, pp. 397-405 (1988). 

In general, the Algebra Reconstruction Technique 
25 requires a large amount of computation time, and has not been 
implemented with parallel processing techniques. 

In view of the parallel processing aspect of the 
present invention, the system described in Richard A. Robb, 
Arnold H. Lent, Barry K. Gilbert, and Aloysius Chu, "The 
30 Dynamic Spatial Reconstructor", J. Med. Syst . , Vol. 4, No. 2, 
pp. 253-288 (1980) is relevant. The Dynamic Spatial Recon- 
structor employs twenty-eight x-ray sources and twenty-eight 
x-ray imaging systems in a synchronous scanning system to 
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acquire data all at one time for a subsequent "stack of 
slices" reconstruction using conventional 2D reconstruction 
algorithms. The Robb et al. reference refers to the use of 
"high-speed parallel processing techniques" for the recon- 
5 struction computation. 

Summary of thg> Tn^*=>ni-4r7p 

Accordingly, it is an object of the invention to 
provide practical methods and apparatus for reconstructing 3D 
images from cone beam projection data, in view of the massive 
computational power required. 
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It is another object of the invention to provide 
methods and apparatus employing parallel processing for 
reconstructing 3D images from cone beam projection data. 

It is another object of the invention to provide 
methods and apparatus which can reconstruct selected regions 
15 of interest/ without spending processor time reconstructing 
regions outside the regions of interest. 

Briefly, and in accordance with an overall aspect 
of the invention, a parallel processing architecture and 
method, based on the 3D Algebra Reconstruction Technique 
(ART) , are defined for iteratively reconstructing from cone 
beam projection data a 3D voxel data set V representing an 
image of an object. The cone beam projection data is 
acquired by employing a cone beam x-ray source and a two- 
dimensional array detector to scan the object along a source 
25 scanning trajectory to obtain, at each of a plurality of 
discrete source positions 0/ on the source scanning tra- 
jectory, a 2D measured cone beam pixel data set , The 
invention is applicable to reconstructing an object image 
from cone beam projection data acquired with any scanning 
30 geometry, so long as it is well defined. Examples include 
source scanning trajectories which comprise a single circle. 
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dual or multiple parallel circles, dual perpendicular cir- 
cles, or square waves on a cylindrical surface. However, in 
the embodiment described in detail hereinbelow, a circular 
source scanning trajectory lying in a plane perpendicular to 
5 a rotation axis passing through the object and perpendicular 
to all planes of the detector is employed. In this case, the 
discrete source positions 9/ are defined simply as angular 

positions. With various other scanning trajectories, the 
discrete source positions G/ are defined by parameters in 
10 addition to or instead of angular position. 

In accordance with the invention, the 3D voxel data 
set V is subdivided into or organized as a plurality of m of 
voxel subcubes through V"^'^ which are independent of each 
other. In other words, there is no overlap between adjacent 

15 voxel subcubes. Each voxel subcube contains a plurality of 
voxels. In the illustrated embodiment, boundaries between 
adjacent voxel subcubes lie in planes which are perpendicular 
to the detector planes and parallel to each other. However, 
independent voxel subcubes may be organized in a variety of 

20 ways. 

As a result of the subdivision of the 3D voxel data 
set V into voxel subcubes, the 2D measured cone beam pixel 
data set P- (measured projection data array) is correspond- 
ingly subdivided for each source position 6/. In the embodi- 

25 ment described in detail herein wherein boundaries between 
adjacent voxel subcubes lie in parallel planes, the measured 
projection data array is divided into what may be termed pro- 
jection data strips. In the more general case, the measured 
projection data array is divided into what may be termed pro- 

30 jection data subsets. 

Thus, for each source position 0/, each voxel sub- 
cube has a corresponding projection data strip or projection 
data subset. However, unlike the voxel subcubes, the projec- 
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tion data strips or subsets are not independent. In other 
words, dependent upon the particular scanning geometry 
defined and the particular voxel and pixel positions for each 
particular source position 0/, there are overlapping regions 
5 between one or more adjacent projection data strips or sub- 
sets corresponding to adjacent voxel subcubes . 

Nevertheless, in accordance with the invention, 
each voxel subcube and its corresponding projection data 
strip or subset can be processed in parallel with other pairs 
of voxel subcubes and corresponding projection data strips or 
subsets, without interference. A bi-level parallel-processor 
architecture is employed. Backpro ject /repro ject processor in 
level 0 selectively perform voxel driven backpro jection from 
the projection data strips or subsets to the corresponding 
15 voxel subcubes, and voxel-driven reprojection from the voxel 
subcubes to the corresponding projection data strips or sub- 
sets. Each backpro ject/repro ject processor operates on data 
which is independent of the data for the other backpro- 
ject /repro ject processors, so the processors can operate in 
20 parallel. At least one split/merge processor in level 1 

splits a projection data set for a particular source position 
9i, into projection data strips or subsets from reprojection 
into a projection data set. During iterative reconstruction 
implementing the Algebra Reconstruction Technique, for each 
25 source position 6/, data are sent back and forth between level 
0 for backprojection and reprojection and level 1 for splitt- 
ing and merging. 

Preferably, level 1 is organized in a tree-type 
manner as a plurality of layers. The split/merge processors 
30 can comprise transputers. 



More particularly, a method in accordance with the 
invention includes the steps of organizing the voxel data set 
V as a plurality m of independent voxel subcubes through 
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V"^'^ f each voxel subcube including a plurality of voxels; ini- 
tializing the voxel data set V; and performing successive 
iterations to correct the values of at least selected voxels 
of the 3D voxel data set based on the measured cone beam data 

5 sets . In a particular embodiment, boundaries between adja- 
cent voxel subcubes lie in planes perpendicular to all planes 
of the 2D detector array and parallel to each other. 

During each iteration, the values of at least the 
selected voxels for each of the discrete source positions 0/ 

10 on the source scanning trajectory are corrected by reproject- 
ing each voxel subcube through V^'^ to calculate pixel val- 
ues for a corresponding 2D calculated projection data strip 
or subset from a group of calculated projection data strips 
or subsets P/^ through Z^"*'^ of a 2D calculated projection data 

15 set Pi including a plurality of pixels, the projection data 

strips of subsets P/* through P/^"^ overlapping in part; merging 
the calculated projection data strips of subsets P/* through 
P^"^ to calculate pixel values of the 2D calculated projection 
data set Pi, the step of merging including adding pixel values 

20 in any overlapping regions of the calculated projection data 
strips of subsets P/^ through P^~^ ; calculating a 2D correction 
projection data set Ei including a plurality of pixels by 
determining the normalized difference between the value of 
each pixel of the measured projection data set P and the 

25 value of the corresponding pixel of the calculated projection 
data set Pi splitting the 2D correction projection data set £/ 

into a plurality of 2D correction projection data strips or 
subsets £f through E*"^ corresponding to the voxel subcubes 
through V^'^ , the 2D correction projection data strips or sub- 
30 sets through E^"^ overlapping in part, including duplicat- 
ing element values for any overlapping regions of the correc- 
tion projection data strips or subsets E^ through Ep"^ ; and 
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backprojecting each correction projection data strip or sub- 
set Ei through E^^^ to correct voxel values in the correspond- 
ing voxel subcube of the plurality of voxel subcubes Vo 
through V^'^ , 

5 The step of merging is organized in a tree-type 

manner as a plurality of layers and includes, for each layer, 
merging groups of calculated projection data strips or sub- 
sets from the next lower layer to present a lesser number of 
calculated projection data subsets to the next higher layer, 
10 the lowest layer merging groups of the calculated projection 
data strips or subsets through Pf^"^ from the step of repro- 

jecting, and the highest layer producing the calculated pro- 
jection data set Pi. 

The step of splitting is also organized in a tree- 
15 type manner as a plurality of layers and includes, for each 
layer, splitting one or more correction projection data 
strips or subsets from the next higher layer to present a 
greater number of correction projection data strips or sub- 
sets to the next lower layer, the highest layer splitting the 
20 correction data set £/, and the lowest layer producing the 

correction projection data strips or subsets E° through E^'^ . 

The method further includes employing a plurality m 
of backproject/reproject processors corresponding to the 
voxel subcubes to perform the steps of reprojecting and back- 
25 projecting, and employing a plurality of split/merge proces- 
sors connected in a tree structure to perform the steps of 
merging and splitting. 

The method additionally includes, during each iter- 
ation, for each of the discrete source positions G/, in con- 
30 junction with the steps of reprojecting each voxel subcube 
through V»*-^^ calculating data values for a corresponding 2D 
normalizing factor data strip or subset from a group of 2D 
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normalizing factor strips or subsets through A^"*^ of a 2D 
normalizing factor data set Ni including a plurality of pix- 
els, the 2D normalizing factor strips or subsets through 
A^*"^ overlapping in part; merging the normalizing factor data 
5 strips or subsets through N^'^ to calculate element values 
of the 2D normalizing factor data set and, during the step 
of calculating the 2D correction data set £/, employing the 

value of each corresponding element of the normalizing factor 
data set iV/. 

10 The step of reprojecting is voxel driven and 

includes, for each voxel subcube of the plurality of voxel 
subcubes through V^'^ , initializing the corresponding cal- 
culated projection data strip or subset from the group of 
projection data subsets P** through initializing the ccr- 

15 responding normalizing factor data strip or subset from the 

group of normalizing factor data strips or subsets through 
N""^; and for each of the selected voxels of the 3D voxel data 

set V included in the particular voxel subcube, adding to the 
value of each pixel of the corresponding calculated projec- 

20 tion data strip or subset from the group of calculated pro- 
jection data strips or subsets through /^"""^ influenced by 
the selected voxel the voxel value multiplied by a weighting 
factor h{p,v) determined for the voxel and pixel positions 
based on geometry such that pixel values are cumulated in the 

25 corresponding artificial projection data subset, and adding 
to the value of each element of the corresponding normalizing 
factor data strip or subset from the group of normalizing 
factor data strips or subsets through N^'^ the square of 

the weighting factor h(p,v) determined for the voxel and pixe 
30 positions based on geometry such that element values are 

cumulated in the corresponding normalizing factor data strip 
or subset- 
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The step of backpro jecting is also voxel driven and 
includes, for each voxel subcube of the plurality of voxel 
subcubes Vo through V»«-^ for each of the selected voxels of 
the 3D voxel data set V included in the particular voxel sub- 
5 cube, adding to the voxel value the value of each element of 
the corresponding 2D correction projection data strip or sub- 
set of the group of correction projection data strips or sub- 
sets through £""^ in a pixel position influenced by the 
selected voxel multiplied by a weighting factor h(p,v) deter- 
10 mined for the voxel and pixel positions based on geometry 
such that corrected voxel values are cumulated. 

Parallel processing apparatus implementing the 
Algebra Reconstruction Technique in accordance with the 
invention includes a data memory for storing the voxel data 

15 set V organized as a plurality m of independent voxel sub- 
cubes VO through V^-^ , each voxel subcube including a plural- 
ity of voxels; and processing elements for initializing the 
voxel data set V and performing successive iterations to cor- 
rect the values of at least selected voxels of the 3D voxel 

20 data set based on the measured cone beam data sets during 

each iteration correcting the values of at least the selected 
voxels for each of the discrete source positions 6/ on the 
source scanning trajectory. 

The processing elements include in a level 0 a plu- 
25 rality m of backpro ject/reprojectors corresponding to the 
voxel subcubes, and in a level 1 at least one split /merge 
processor. The backpro ject/reproject processors are operable 
to reproject each voxel subcube Vo through V^'^ to calculate 
pixel values for a corresponding 2D calculated projection 
30 data subset from a group of calculated projection data sub- 
sets P/" through P/^'^ of a 2D calculated projection data set Pi 
including a plurality of pixels, the projection data subsets 
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through P^^^ overlapping in part. The split/merge proces- 
sor is operable to merge the calculated projection data sub- 
sets Pi through P^'^ to calculate pixel values of the 2D 
calculated projection data set Pi, adding pixel ^values in any 
5 overlapping regions of the calculated projection data subsets 
Pi through P^"^ • The processing elements include means for 
calculating a 2D correction projection data set £i including 
a plurality of pixels determining the normalized difference 
between the value of each pixel of the measured projection 
10 data set pixel P^ and the value of the corresponding pixel of 
the calculated projection data set Pi, The split/merge pro- 
cessor is also operable to split the 2D correction projection 
data set £/ into a plurality of 2D correction projection data 

subsets El through E""^ corresponding to the voxel subcubes V 
15 through V^'^ , duplicating element values for any overlapping 
regions of the correction projection data subsets El through 
iE^""^ - The backpro ject/repro ject processors are also operable 
to backproject each correction projection data subset El 
through E^""^ to correct voxel values in the corresponding 
20 voxel subcube of the plurality of voxel subcubes through 

Preferably there are a plurality of split/merge 
processors organized as a tree-type structure in a plurality 
of layers. The split/merge processors are operable to merge 
25 groups of calculated projection data subsets from the next 

lower layer to present a lesser number of calculated projec- 
tion data subsets to the next higher layer, the lowest layer 
merging groups of the calculated projection data subsets P^ 
through fj"""^ from the backpro ject/repro ject processors , and 

30 the highest layer producing the calculated projection data 

set Pi. The split/merge processors are also opercUDle to split 

one or more correction projection data subsets from the next 
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higher layer to present a greater number of correction 
projection data subsets to the next lower layer, the highest 
layer splitting the correction data set £/, and the lowest 
layer producing the correction projection data subsets E- 
5 through £;""^ 

At least the split/merge processors comprise trans- 
puters. The split/merge processors of each layer each 
include on I/O connected to the next higher layer, and three 
I/O parts connected to the next lower layer. 
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Brief Description nf i-h^ Draw-inrjrcs 



While the novel features of the invention are set 
forth with particularity in the appended claims, the inven- 
tion, both as to organization and content, will be better 
understood and appreciated, along with other objects and fea- 
tures thereof, from the following detailed description taken 
15 in conjunction with the drawings, in which: 

FIG. 1 depicts data acquisition employing cone beam 
geometry for scanning an object to acquire cone beam projec- 
tion data, and additionally graphically depicts the geometry 
for backprojection and reprojection operations; 

2° FIG. 2 represents the steps in the prior art 3D 

Algebra Reconstruction Technique (ART) for reconstruction; 

FIG. 3 depicts the organization of a 3D voxel data 
set into independent voxel subcubes; 

FIG. 4 depicts the geometry relationship between 
25 voxel subcubes and projection strips or subsets; 



FIG. 5 is a black diagram of apparatus implementing 
ART in accordance with the invention; 
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FIG. 6 represents the steps of reprojection and 
merge in the ART implementation of the invention; 

FIG. 7 represents the steps of split and backpro- 
ject in the ART implementation of the invention; 

5 FIG. 8 depicts the merging of projection strips or 

subsets; and 

FIG. 9 depicts splitting into projection strips or 

subsets . 

Detailed Description 

Referring first to FIG. 1, depicted in a typical 
10 scanning and data acquisition configuration employing cone 
beam geometry. A cube 20 alternatively represents an actual 
3D object being scanned or a 3D voxel (volume element) data 
set V representing the actual object. The voxel data set V 
may be indexed, for example, by voxel number v. It will be 
15 appreciated that, during a scanning and data acquisition 
operation, the actual object is scanned with x-rays; and, 
during a subsequent image reconstruction operation, the voxel 
data set V is calculated to represent an image of the actual 
object. The present invention is directed to image recon- 
20 struction. 

A global coordinate system Cy=(X^, F^„ZJ is defined, 
with its origin (0,0,0) at the center of the voxel data set V 
or object 20. The object 20 is positioned within a field of 
view between a cone beam x-ray point source 22 and a 2D array 
25 detector 24, from which cone beam projection data is 

acquired. The 2D array detector 24 is connected to a conven- 
tional high speed data acquisition system (not shown) . An 
axis of rotation 26 coincident with the Zy axis passes through 

the field of view and the object 20. Perpendicular to the Zy 
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axis of rotation 2 6 is a midplane within which the Xy and Yy 
axes lie, as well as a circular source scanning trajectory 28 
centered on the Zy axis of rotation 26. 

For scanning the object at a plurality of discrete 
source position Qi on the source scanning trajectory 28, 
which in FIG. 1 are discrete angular positions 9,-, the x-ray 
source 22 moves relative to the object 20 along the source 
scanning trajectory 28, while the 2D array detector 24 
remains fixed in position relative to the source 22. Either 
the source 22 and detector 24 may be rotated around a sta- 
tionary object 20, or the object 20 may be rotated while the 
source 22 and detector 24 remain stationary. 

At each discrete source position S/ a 2D measured 
cone beam pixel data set ^. is acquired by the data acquis i- 
15 tion system (not shown) connected to the array detector 24, 

and stored for subsequent reconstruction. Each measured cone 
beam pixel data set includes a plurality of pixels indexed, 
for example, by pixel number p. 
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20 



Thus, in the particular geometry of FIG. 1, data 
are acquired at a number of angular positions or view angles 
9/ around the object. As depicted, a is the source-to-detec- 
tor distance and P is the source-to-rotation-axis distance. 
For each view angle 9j, the x-ray source 22 is located at 

[^sin0,,-^cosa„Of and a 2D projection ^ is acquired. The cen- 
25 ter of the 2D projection array 24 is located at 
0, = [-(a-)5)sin0,. (a-j8)cosd,.of . 

FIG. 2 summarizes the procedure of the prior art 
three-dimensional Algebra Reconstruction Technique (ART) , 
which is an iterative procedure using raysums to correct the 
30 value of each voxel of the 3D voxel data set V at each itera- 
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tion. After a number of iterations, the values of the voxel 
data set V converge to represent the actual 3D object. 

In box 30, the procedure commences by initializing 
the 3D voxel data set V, for example by clearing all voxel 
5 values to zero. The procedure then goes through successive 
iterations and, for each iteration, through each source posi- 
tion or projection angle 9i corresponding to discrete angular 

positions on the source scanning trajectory 28. For each 
projection angle 9/ there are three steps: repro jection, 

10 error calculation and backpro jection. The procedure is com- 
pleted (not shown) either when a predetermined number of 
iterations have been processed, or when the error is below a 
predetermined threshold. The projection angles 6/ may be 

processed in any order, either sequentially or in some other 
15 order Depending upon the particular scanning geometry, the 

order of projection angle 9/ processing can affect the rate of 

convergence of the voxel data set V. 

Box 32 states the step of repro jection, which is in 
turn detailed in Box 34 in the manner of a subroutine. The 
20 3D voxel data set V (or a selected portion of interest) is 

reprojected from real space to projection space to generate a 
2D calculated projection data set Pi for the particular pro- 
jection angle 9|. The 2D calculated projection data set 

corresponds on a pixel-by-pixel basis with the actual mea- 
25 sured cone beam pixel data set P^, but with actual pixel data 

values differing depending upon how accurately the voxel data 
set V represents the actual object, and depending upon vari- 
ous inherent computational errors. In conjunction with the 
repro jection, a 2D normalizing factor data set Ni is also 

30 calculated, corresponding on a pixel-by-pixel basis with the 
calculated and measured projection data sets P/ and . 

Considering the reprojection process of Box 34 in 
detail, the calculated projection data set Pi and the normal- 



10 
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izing factor data set iV/ are each initialized by clearing all 
pixels to zero. Then for each voxel v of the 3D voxel data 
set V, each pixel p of the calculated projection data set Pi 
and the corresponding normalizing factor data set iV,- influ- 
5 enced by the voxel is identified in a manner described here- 
inbelow with reference again to FIG. 1. For each pixel p so 
influenced by voxel v, the following three operations are 
performed: (1) A weighting factor h(p,v) for the particular 
voxel and pixel positions is determined based on geometry, 
for example by a well-known bi-linear interpolation process. 
(2) The value of the voxel V(v) is multiplied by the weight- 
ing factor h(p,v) and added to the value of the calculated 
projection data set pixel Pi(p) such that pixel values are 
cumulated. (3) The weighting factor h(p,v) is multiplied by 
15 itself and added to the value of the normalizing factor dat 
set element N^p) such that element values are cumulated. 

Returning to Box 30 of FIG. 2, a 2D correction pro- 
jection data set £i indexed, for example, by pixel number p, 
is calculated. For each pixel of the corresponding data 
sets, the difference between the value of the measured pro- 
jection data set pixel Pi(p) and the value of the corresponding 
calculated projection data set pixel Pi(p) is determined, and 
then normalized by dividing by the value of the corresponding 
normalizing factor data set element Ni(p) . The result is the 

25 value of the corresponding correction projection data set 
element £ifp>. 



20 



Box 36 states the step of backpro jection, which is 
in turn detailed in Box 38 in the manner of a subroutine. 
The 2D correction projection data set for the particular 
projection angle 9,- is backpro jected from projection space to 

real space to correct voxel values in the 3D voxel data set 
V. 
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Considering the backpro jection process of Box 38 in 
detail, for each voxel v of the 3D voxel data set V, each 
pixel p of the correction projection data set £/ influenced by 

the voxel is identified. For each pixel p so influenced by 
5 voxel V, the weighting factor h{p,v) is determined via inter- 
polation based on geometry, and then the value of the correc- 
tion projection data set pixel Ei(p) is multiplied by the 

weighting factor h(p,v), and the product added to the value 
of the voxel value V(v), such that voxel values are cumula- 
10 tively corrected. 

Referring again to FIG. 1, the manner in which 
array detector pixels influenced by each voxel is determined 
will now be described in the context of the backpro jection 
and repro jection procedures. 

15 Let V(Xy, y^, Zy) represent a three-dimensional voxel 

data set defined in the three-dimensional discrete Cartesian 
coordinate system, C^, with the origin located at the center 
of the 3D cubic data. For a 3D voxel data set with 
{Nx^NyXNz) voxels, every voxel can be indexed by (i,,/^,,/,). 

20 The coordinate of each voxel is stated as 

where Kx, Ky, and Kz denote the voxel resolution along the x, 
y, and z axis in the coordinate system Cy, respectively. A 

two-dimensional discrete Cartesian coordinate system C^^ is 

25 used to address each detector element in the 2D area detec- 
tor. Every detector element D{xci/ yd) is indexed by 
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(2) 



the tlx and \ly are used to indicate the detector resolution 
along the x and y axis, respectively. The relative relation- 
ship between Cp^ and Cl can be stated as 4x4 transformation 
function Tp^_^. For each detector Pa.iXj,yj) defined in the 
detector coordinate system Cj^, the coordinate (Xd, yd) can be 

transferred to the voxel coordinate system by applying the 
transformation function 7>^_v,i.e., 



yd 

0 



(3) 



10 



transformation function Tp^^,, which transfers the detector 

coordinate system to the voxel coordinate system, is 
expressed as rotating about the y axis by 6/, then rotating 

90* around the reference x axis, and finally translating by 
Op; i.e.. 
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(4) 



where 0/ is View angle measured respect to the y axis in Cl 

a is Source to detector distance (STOD) 

P is Source to rotation axis distance (STOD) 

then Ecjuation (3) can be rewritten as 

'cos d 0 sin ©i - (a - /?)sin 



sin©. 0 -cos©,(a-i3)cos©i 
0 10 0 
0 0 0 1 



0 



(5) 



1 .e, 



X = cos ©, - (a - )8)sin ©, 
y = sin 6i + (a - j3)cos 



(6) 



For each voxel V(Xy;, yy, Zy) , given a view angle Of, an integra- 
10 tion line such as an exemplary integration line 44 passes 
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through the x-ray source 22 located at [^sin0., -j3 cos0.,of and 
voxel V(xv,yy,Zy) , can be expressed as 

x-p^md, _ y-bficosd, z 

where (i,y,z) is any point belonging to the integration line 
5 44. This integration line 44 intersects the detector plane 
at the location (x^.y^) in the projection coordinate system 

Cp^. Substituting Equation (6) into Equation (7) yields 

x^ cos - (g - /?) sin g, - p sin g, ^ x^ sin 6, + (a - /?) cos H- cos ft 
jc, ^ jSsinft +jScosft " 

Then the intersection point (xd/ yd) on the detector or pro- 
10 jection plane can be calculated as 

^ _ Q:[x,cosg,-Hj^^sing,1 
y^ cos ft - X, sin ft. + j3 

az 

y^ 2s ^ 

j'vCosft.-jc^sinft+ZJ 

A voxel-driven approach is applied for both back- 
projection and reprojection. At a given view angle Qi, a 

15 unique intersection point (Xd, yd) on the detector plane C^^ 
can be calculated for a given voxel (Xvr yv/ Zv) in , as 
represented in FIG. 1 by the intersection point 40 and voxel 
42 along line 44. Four immediate neighbor pixels surrounded 
the intersection point (xdr yd) • By a known process of bi- 
linear interpolation, weighting factors reflecting the con- 
tribution of each voxel to each pixel on the detector plane 
can be calculated for reprojection, and weighting factors 
reflecting the contribution of each pixel on the detector 
plane to each voxel can be calculated for backpro jection. 



20 
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With reference now to FIG. 3, the 3D voxel data set 
V (or 20 in FIG. 1) is organized as a plurality m of indepen- 
dent voxel subcubes through V^'^ , The voxel subcubes are 
independent in the sense that no overlapping re'gion exists 
5 between adjacent voxel subcubes. While the voxel data set 
may be divided in a variety of ways, preferably boundaries 
between adjacent voxel subcubes lie in planes which are per- 
pendicular to the detector planes and parallel to each other, 
as depicted in FIG. 3. The size of each voxel subcube is 
10 constrained by the memory space associated with each proces- 
sor element used to do back projection and repro jection, 
described hereinbelow with reference to FIGS. 5, 8 and 9. 

FIG. 4 depicts the geometric relationship between 

the voxel subcubes of FIG. 3 and projections on the detector 
15 array for each angular position S/. For each voxel subcube in 

real space there is a corresponding projection subset on the 
projection plane. In the planar arrangement of voxel sub- 
cubes illustrated, the projection subsets take the form of 
parallel strips on the detector plane. 

20 Thus, in FIG. 4 there are two representative voxel 

subcubes 46 and 48, with corresponding projection strips 50 
and 52. There is an overlapping region 54 between the two 
projection strips 50 and 52. It is however a feature of the 
invention that each voxel subcube and its corresponding pro- 

25 jection strip can be projected and reprojected in parallel 
with other pairs of voxel subcubes and projection strips 
without interference . 

FIG. 5 depicts apparatus in accordance with the 
invention, comprising a bi-level parallel processor architec- 
30 ture. Included in the FIG, 5 apparatus is a data memory 60 
for storing the 3D voxel data set V, organized as indicated 
into separate memory areas for each of the independent voxel 
subcubes through V^'^ . 
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In a level 0 are a plurality of m of backpro- 
ject/repro ject processors, such as processor 62, correspond- 
ing respectively to the voxel subcubes , Each of the backpro- 
ject/reproject processors of level 0 operates independently 
5 of an in parallel with the other processors of level 0 to 

reproject and backproject back and forth between a particular 
voxel subcube and a particular projection strip as part of 
the overall iterative process. 

In a level 1 is at least one, and preferably many, 
10 split/merge processor, such as processor 64, Using 4-link 
connectable processors, such as the INMOS T800 transputer, 
processors in level 1 are organized in layers in a tree-type 
manner. Each level 1 processor (or node) has three sons and 
one parent in the particular embodiment illustrated. The 
15 INMOS T800 transputer is basically a computer on a chip with 
a 32-bit integer processor and a 64 bit floating point pro- 
cessor,r with a practically achievable floating point computa- 
tion rate of approximately 1.5 MFLOPS. It includes 4 kBytes 
of on-chip memory, and address space for 4 Gbytes of off-chip 
20 memory. Each transputer can be connected to four neighbors 
over a point-to-point link with a 20-Mbit /second bandwidth. 
An available development system, namely a Meiko Computing 
Surface, includes 512 transputers, with a maximum computation 
rate of approximately 800 MFLOPS. 

25 The level 1 processors are operable, following a 

level 0 reproject operation, which produces m calculated pro- 
jection data subsets through to merge the calculated 
projection data subsets through /J-'""^ into a memory 66, as 
well as to merge normalizing factor data subsets iVf through 

30 iV/""^ from the level 0 processors into a single normalizing 

factor data set Ni for storing in the memory 66. The level 1 
processors are also operable to split a correction projection 
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data set £/ stored in the memory 66 into a plurality of cor- 
rection projection data subsets through for backpro- 
jection by the level 0 processors. In FIG. 5, the arrows 68 
and 70 respectively indicate these two types of operation, 
5 namely, reproject then merge, and split then backproject. 

Also shown in FIG. 5 is a correction processor 72 
which operates in a conventional manner as described herein- 
above to calculate the correction projection data set £"/ based 
on the difference between the calculated projection data set 
10 Pi, and the measured projection data set ^ for the particular 
projection angle 0i- The measured projection data sets ^- 
resuiting from data acquisition as described hereinabove with 
reference to FIG. 1 are stored in exemplary memories 74, 7 6, 
78 and 80. 

15 Thus,, one projection array is associated with each 

split/merge processor in level 1. During backpro jection, 
each split/merge processor reads in a projection reads in a 
projection strip from its parent node, splits it into three 
projection substrips, and then distributes those to its son 

20 nodes for further splitting or backpro jection . During repro- 
jection, each split/merge processor reads in three projection 
strips from its son nodes (either the backpro ject/reproject 
processors or split/merger processors) , merges them into one 
projection strip, and sends that projection strip to its par- 

25 ent node for further merging. The number of processor layers 
in level 1 is related to the total number of backpro- 
ject/reproject processors in level 0. One voxel subcube and 
its corresponding projection strip are associated with each 
backpro ject/reproject processor in level 0. During backpro- 

30 jection, each split/merge processor reads in a projection 
strip from its parent node and splits it into three projec- 
tion substrips, and then distributes those to its son nodes. 
Each backproject/reproject processor receives a projection 
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strip from its parent node and performs voxel-driven backpro- 
jection. During repro jection, each backpro ject/repro ject 
processor calculates a projection strip for the desired view 
angle and sends the resulting projection strip to its parent 
5 node (split /merge process -r) for merging. The 3D algebra 

Reconstruction Technique is achieved by iteratively perform- 
ing backpro jection and repro jection. The entire parallel pro- 
cessing procedure of the 3D cone beam algebra reconstruction 
algorithm can be divided into four major tasks: split, 
10 merge, backpro ject, and repro ject. 

These four tasks are represented in Boxes 82 and 84 
of FIGS . 6 and 7 . In the ART implementation of the inven- 
tion, the repro jection and merge tasks in Box 82 of FIG. 6 
are substituted for the prior art reprojection step repre- 
15 sented in Box 34 of FIG. 2. Similarly, the split and back- 
project tasks in Box 84 of FIG. 7 are substituted for the 
prior art backpro jection step represented in Box 38 of FIG. 
2. The four tasks will now be considered in detail in the 
order presented in the Boxes 82 and 84. 

20 The reproject process is executed by each backpro- 

ject/repro ject processor in level 0 or reprojection. Its 
primary function is to generate a projection strip for a 
given view angle from a voxel subcube and send the results to 
its parent node for projection merging. The reprojection 

25 procedure uses a voxel-driven approach and is illustrated as 
follows : 
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for a given view angle 6i 

set P = 0 and N = 0 

for each voxel v in the voxel subcube 

calculate pixel location in the projection strip influenced by v for each 
5 pixel p influenced by v 

calculate h(v,p), the weighting factor, via interpolation 
P(p) = P(p) + h(p,v)V(v) 
N(p) = N(p) + h(p.v)h(p.v) 

end for 

10 endfor 
end for 

The merge process runs on the split /merge proces- 
sors in level 1 for reprojection purposes. Its primary func- 
tion is to merge three projection strips, read from its three 

15 son nodes, into one. As mentioned hereinabove, an overlap- 
ping region between two projection strips corresponds to two 
adjacent voxel subcubes . The merge is achieved by adding 
projection strips together in the overlapping region accord- 
ing to its global coordinate system, as represented in FIG. 

20 8. 

The split process is executed by the split/merge 
processors in level 1 for backpro jection . Its primary func- 
tion is to partition the projection strip into three projec- 
tion substrips according to the global location of the voxel 

25 subcubes processed by its son nodes. As mentioned herein- 
above, an overlapping region exists between projection strips 
that corresponds to two adjacent voxel subcubes. Projection 
data in the overlapping region are duplicated in two adjacent 
projection strips, as represented in FIG. 9. Before perform- 

30 ing projection splitting, the processor reads in projection 
Strips from its parent node and reads in information from its 
son nodes about the starting point and dimensions of the 
voxel sector of its son nodes. Then the input projection 
strip is divided into three projection substrips, according 

35 to the global location of the voxel subcubes processed by its 
son nodes and the scanning geometry, and the resulting three 
projection strips are distributed to its son nodes. 
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The backproject process is executed by the backpro- 
ject/reproject processors in level 0 for backpro jection . Its 
primary function is to read in a projection strip for every 
view angle and backproject it to a voxel subcube. The back- 
5 projection procedure uses a voxel-driven approach and is 
illustrated as follows: 

for every view angle Gj 

lead in projection strip from its parent node 
for each voxel v in the voxel subcube 
^ ^ calculate pixel location in the projection strip influenced by v for each 

pixel p influenced by v 

calculate h(v,p), the weighting factor, via interpolation 
V(v) = V(v) + h(v.p)P(v) 
- c N(p) = N(p) + h(p,v)h(p.v) 

endfor 

end for 

endfor 

Considering the parallel processing approach in 
detail with reference to the geometry of FIGS. 1, 3 and 4, an 
(N^xNyXN^) voxel data set can be divided into (M) voxel sub- 
cubes, each with [({N^/M)xN;^xN^)] voxels. The starting pint 
of each voxel subcube is stated as 



20 



<:y]J A^z r^zl 1 



for m=0, . . (M-1) . The variable m is used to denote the 
25 voxel subcube index. For each voxel subcube, its correspond- 
ing region on the projection plane that may contribute to it 
can be identified according to the scanning geometry. There 
are M projection strips corresponding to voxel subcubes . The 
starting point of the projection strip corresponding to the 
30 mth voxel subcube is stated as 
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and each projection strip contains 



[^ + «/2)] 



projection elements. The variable denotes the number of 



detector elements along the x axis of the projection coordi- 
5 nate system. The geometric relationship between the starting 
point of each voxel subcube and the starting point of its 
corresponding projection strip is illustrated in FIG, 4, 

Each voxel in a voxel subcube handled by a backpro- 
ject/reproject processor is indexed by {iv»jv»/^v)f where 



The voxel resolutions along the x, y, and z axes are speci- 
fied as {Ry.R^^Ry). The integration line connecting the x-ray 
15 source and the voxel (Vjf,Vy,V^) for the view angle 9i can be 
stated as 



This integrating line will intersect the detector plane at 
the location (X^, Yd) in the projection coordinate system show 

20 in FIG. 4. 



a = o,...,(n;-i) 



10 




Its global coordinate can be calculated as 
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^ _ «[v;rCOsg.. + yySing,l 
V;.cosa.-V<fSin0i+^ 



In order to reduce computation, the {Xd,Yd) can be calculated 
as follows : 

X, = -?^ 
" B 
y _ aC 

where 



^ = A) + (^^ cos 0. ) + y,, (/?^ sin 0, ) 
^ = ^0 - sin d,) + cos a,) 



and 



A = SlR^ cos d, + S^;?^ sin 6., 
B = S^R} cos d, - 5;i[/?^f sin 0, + ^ 



Cq — S2RV 



10 



15 



20 



The values A, B, and C can be calculated by incrementing 
their values with a constant increment. This approach sig- 
nificantly reduces the computation complexity of calculating 
Xd and Yd. Then a bi-linear interpolating scheme is applied 
to calculate the weighting factor for all the pixels influ- 
enced by the voxel {Vx. Vy, Vz) . The computation of the pixel 
location (Xd.Yd) is unique and applicable to both the back 
projecting and the reprojection cases. 

The bi-level tree-type parallel-proessor architec- 
ture and the four major processing tasks (split, merge, back- 
project, and reproject) are implemented on the Meiko Comput- 
ing Surface (T800 transputer development system) . Each back- 
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pro ject/repro ject processor of each split/merge processor 
(FIG. 5) is mapped to a T800 transputer. 

While the illustrated embodiment employs a geometry 
wherein there is overlap only between .two adjacent projection 
5 strips or subsets, it will be appreciated that more than two 
projection strips or subsets may overlap. Similarly, each 
FIG, 5 split/merge processor in level 1 has one parent node 
and three son nodes, which facilitates implementation with 
transputers which have four I/O parts. However, different 

10 numbers of son nodes may be employed. Moreover, rather than 
a three-type architecture^ level 1 may comprise a single pro- 
cessor which merges all projection strips or subsets in a 
single operation, or correspondingly splits the complete pro- 
jection data set into a plurality of projection strips or 

15 subsets in one operation. 

This technique can be applied to both medical and 
industrial 3D CT imaging. Because it is based on a voxel- 
driven algebraic reconstruction technique, it is particularly 
useful for reconstructing regions of interest. The paral- 
20 lelism provided by the Split and Merge Reconstruction method 
can be extended to large data sets . Performance improves as 
the number of backpro ject/repro ject processors grows. 

While specific embodiments of the invention have 
been illustrated and described herein, it is realized that 
. 25 modifications and changes will occur to hose skilled in the 
art . It is therefore to be understood that the appended 
claims are intended to cover all such modifications and 
changes as fall within the true spirit and scope of the 
invention . 
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What is claimed is: 

1. A parallel processing method based on the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 
5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
source scanning trajectory to obtain, at each of a plurality 
of discrete source positions Q± on the source scanning tra- 
jectory, a 2D measured ccr.e beam pixel data set P^, said 
10 method comprising: 

organizing the voxel data set V as a plurality m of 
independent voxel subcubes Vo through Vm-1 , each voxel subcube 
including a plurality of voxels; 

initializing the voxel data set V; and 
performing successive iterations to" correct the 
values of at least selected voxels of the 3D voxel data set 
based on the measured cone beam data sets Pi, during each 
iteration correcting the values of at least the selected 
voxels for each of the discrete source positions 6/ on the 
20 source scanning trajectory by 

reprojecting each voxel subcxibe V through V'"--' 
to calculate pixel values for a corresponding 2D 
calculated projection data subset from a group of 
calculated projection data subsets Pf through P"~^ 
of ^ 2D calculated projection data set Pi including 
a plurality of pixels, the projection data subsets 
P' through P^~^ overlapping in part 

merging the calculated projection data subsets 
through /J.""' to calculate pixel values of the 2D 
calculated projection data set Pi, the step of 
merging including adding pixel values in any over- 
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lapping regions of the calculated projection data 
subsets through Pr'\ 

calculating a 2D correction projection data 
35 set Ei including a plurality of pixels by determin- 

ing the normalized difference between the value of 
each pixel of the measured projection data set ^ 

and the value of the corresponding pixel of the 
calculated projection data set Pi, 

40 splitting the 2D correction projection data 

set Ei into a plurality of 2D correction projection 

data subsets through £7*"^ corresponding to the 

voxel subcubes through v"^-l, the 2D correction 

projection data subsets E? through E"^"^ overlapping 

45 in part, and the step of splitting including dupli- 

cating element values for any overlapping regions 
of the correction projection data subsets 
through E^"^ , and 

backprojecting each correction projection data 

50 subsets through to correct voxel values in 

the corresponding voxel subcube of the plurality of 
voxel subcubes through V^'^ . 

2, A method in accordance with Claim 1, wherein: 

The step of merging is organized in a tree-type 
manner as a plurality of layers and comprises/ for each 
layer, merging groups of calculated projection data subsets 
5 from the next lower layer to present a lesser number of cal- 
culated projection data subsets to the next higher layer, the 
lowest layer merging groups of the calculated projection data 
subsets P^ through P^'^ from the step of reprojecting, and the 
highest layer producing the calculated projection data set Pi; 
10 and 
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the step of splitting is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, splitting one or more correction projection data sub- 
sets from the next higher layer to present a greater number 
15 of correction projection data subsets to the next lower 

layer, the highest layer splitting the correction data set £/, 
and the lowest layer producing the correction projection data 
subsets through E^^ . 

3. A method in accordance with Claim 1, which com- 
prises employing a plurality m of backpro ject/repro ject pro- 
cessors corresponding to the voxel subcubes to perform the 
steps of reprojecting and backpro jecting . 

4. A method in accordance with Claim 2, which com- 
prises employing a plurality of split/merge processors con- 
nected in a tree structure to perform the steps of merging 
and splitting. 

5. A method in accordance with Claim 1, which fur- 
ther comprises, during each iteration, for each of the dis- 
crete source positions 9/: 

in conjunction with the step of reprojecting each 
5 voxel subcube Vo through V^'^ , calculating data values for a 
corresponding 2D normalizing factor data subset from a group 
of 2D normalizing factor subsets through N^'^ of a 2D nor- 
malizing factor data set iV/ including a plurality of pixels, 
the 2D normalizing factor subsets N? through iV""^ overlapping 
10 in part; 

merging the normalizing factor data subsets 

through N^'^ to calculate element values of the 2D normalizing 
factor data set AT/; and 

during the step of calculating the 2D correction 
15 data set £/, employing each corresponding element of the nor- 
malizing factor data set element A^,-. 
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6. A method in accordance with Claim 1, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes 
through V^'^ : 

5 initializing the corresponding calculated projec- 

tion data subset from the group of projection data subsets 

through P^'^ ; and 

for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 
Xo adding to the value of each pixel of the cor- 

responding calculated projection data subset from 
the group of calculated projection data subsets P," 

through P/^'^ influenced by the selected voxel the 
voxel value multiplied by a weighting factor h(p,v) 
15 determined for the voxel and pixel positions based 

on geometry such that pixel values are cumulated in 
the corresponding artificial projection data sub- 
set • 

7. A method in accordance with Claim 5, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes 
through V^'^ : 

5 initializing the corresponding calculated projec- 

tion data subset from the group of projection data subsets P^ 
through P/^'^ ^ 

initializing the corresponding normalizing factor 
data subset from the group of normalizing factor data subsets 
10 Ni through and 

for each of the selected voxels of the 3D voxel 
data set V{v) included in the particular voxel subcube, 

adding to the value of each pixel of the cor- 
responding calculated projection data subset from 
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°f calculated projection data subsets 
through /^"- influenced by the selected voxel and 
voxel value multiplied by a weighting factor h{p,v) 
determined for the voxel and pixel positions based 
on geometry such that pixel values are cumulated in 
the corresponding artificial projection data sub- 
set, and 

adding to the value of each element of the 
corresponding normalizing factor data subset from 
the group of normalizing factor data subsets N; 
through Nr' the square of the weighting factor 
h(p,v) determined for the voxel and pixel positions 
based on geometry such that eLem.«tvt values axe 
cumulated in the corresponding normalizing factor 
data subset. 

8. A method in accordance with Claim 1, wherein 
the step of backprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes Vo 
through V*"'^, 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data subset of the group of correction projec- 
tion data subsets through in a pixel posi- 
tion influenced by the selected voxel multiplied by 
a weighting factor h{p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

9. A method in accordance with Claim 6, wherein 
the step of backprojecting is voxel drivel and comprises, for 
each voxel subcube of the plurality of voxel subcubes Vo 
through VW--?, 
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5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data subset of the group of correction projec- 
10 tion data subsets through £;""' in a pixel posi- 

tion influenced by the selected voxel multiplied by 
the weighting factor h(p,v) determined for the 
voxel and pixel positions based on geometry such 
that corrected voxel values are cumulated. 

10. A parallel processing method based on the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 
5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
circular source scanning trajectory lying in a plane perpen- 
dicular to a rotation axis passing through the object and 
perpendicular to all planes of the detector to obtain, at 
each of a plurality of discrete angular positions 6/ on the 
source scanning trajectory, a 2D measured cone beam pixel 

A 

data set P^, said method comprising: 

organizing the voxel data set V as a plurality m of 
independent voxel subcubes Vo through V^-l , boundaries between 
15 adjacent voxel subcubes lying in planes perpendicular to all 
planes of the detector and parallel to each other, and each 
voxel subcube including a plurality of voxels; 

initializing the voxel data set V; and 
performing successive iterations to correct the 
values of at least selected voxels of the 3D voxel data set 
based on the measured cone beam data sets ^, during each 
iteration correcting the values of at least the selected vox- 
els for each of the discrete angular positions 6/ on the 
source scanning trajectory by 
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25 reprojecting each voxel subcube through V^'^ 

to calculate pixel values for a corresponding 2D 
calculated projection data strip from a group of 
calculated projection data strips through P/^'^ of 
a 2D calculated projection data set Pi, the projec- 

30 tion data strips P/" through P/^'^ overlapping in 

part, 

merging the calculated projection data strips 
P^ through P^'^ to calculate pixel values of the 2D 
calculated projection data set Pi, the step of 
35 merging including adding pixel values in any over- 

lapping regions of the calculated projection data 
strips P^ through , 

calculating a 2D correction projection data 
set Ei by determining the normalized difference 
^9 between the value of each pixel of the measured 

projection data set and the value of the corre- 
sponding pixel of the calculated projection data 
set Pi, 

splitting the 2D correction projection data 
set Ei into a plurality of 2D correction projection 
data strips E^ through EP"^ corresponding to the 
voxel subcubes Vo through V^'^ , the 2D correction 
projection data strips £/ through E^^^ overlapping 
in part, and the step of splitting including dupli- 
cating element values for any overlapping regions 
of the correction projection data strips through 
£^"^ and 

backprojecting each correction projection data 
strip E? through Ep'^ to correct voxel values in the 

55 corresponding voxel subcube of the plurality of 

voxel subcubes through V^'K 
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11. A method in accordance with Claim 10, wherein: 

the step of merging is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, merging groups of calculated projection data strips 
5 from the next lower layer to present a lesser number of cal- 
culated projection data strips to the next higher layer, the 
lowest layer merging groups of the calculated projection data 
strips Z^** through P^'^ from the step of repro ject ing, and the 
highest layer producing the calculated projection data set Pi; 
10 and 

the step of splitting is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, splitting one or more correction projection data- 
strips from the next higher layer to present a greater number 
15 of correction projection data strips to the next lower layer, 
the highest layer splitting the correction data set £/, and 

the lowest layer producing the correction projection data 
strips E- through £/""^ 

12. A method in accordance with Claim 10, which 
comprises employing a plurality m of backpro ject/repro ject 
processors corresponding to the voxel subcubes to perform the 
steps of repro jecting and backpro jecting. 

13. A method in accordance with Claim 11, which 
comprises employing a plurality of split/merge processors 
connected in a tree structure to perform the steps of merging 
and splitting. 

14. A method in accordance with Claim 10, which 
further comprises, during each iteration, for each of the 
discrete source positions 9/; 

in conjunction with the step of reprojecting each 
5 voxel subcube through V^"^, calculating data values for a 
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corresponding 2D normalizing factor data strip from a group 
of 2D normalizing factor strips Nf through N^' of a 2D nor- 
malizing factor data set iV/ including a plurality of pixels, 
the 2D normalizing factor strips through Np"' overlapping 
10 in part; 

merging the normalizing factor data strips 

through NP'^ to calculate element values of the 2D normalizing 
factor data set and 

during the step of calculating the 2D correction 
15 data set £/, employing each corresponding element of the nor- 
malizing factor data set iV/. 

15. A method in accordance with Claim 10, wherein 
the step of rep rejecting is voxel drivel and comprises, for 
each voxel subcube of the plurality of voxel subcubes V» 
through V^'^ : 

5 initializing the corresponding calculated projec- 

tion data strip from the group of projection data strips 
through P/^'^ ; and 

for each of the selected voxels of the 3D voxel 
data set V included in the particular voxel subcube, 

. adding to the value of each pixel of the cor- 
responding calculated projection data strip from 
the group of calculated projection data strips Pp 
through /^""^ influenced by the selected voxel and 
voxel value multiplied by a weighting factor h(p,v) 
determined for the voxel and pixel positions based 
on geometry such that pixel values are cumulated in 
the corresponding artificial projection data strip. 

16. A method in accordance with Claim 14, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes 
through V^-^: 
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5 initializing the corresponding calculated projec- 

tion data strip from the group of projection data strips 
through P^^^ ; 

initializing the corresponding normalizing factor 
data strip from the group of normalizing factor data strips 
10 through N^"^ ; and 

for each of the selected voxels of the 3D voxel 
data set V{v) included in the particular voxel subcube, 

adding to the value of each pixel of the cor- 
responding calculated projection data strip from 

15 the group of calculated projection data strips P^^ 

through P^'^ influenced by the selected voxel and 
voxel value multiplied by a weighting factor h(p,v) 
determined for the voxel and pixel positions based 
on geometry such that pixel values are cumulated in 

20 the corresponding artificial projection data strip, 

and 

adding to the value of each element of the 
corresponding normalizing factor data strip from 
the group of normalizing factor data strips A/,*" 
25 through N^'^ the square of the weighting factor 

h(p,v) determined for the voxel and pixel positions 
based on geometry such that element values are 
cumulated in the corresponding normalizing factor 
data strip. 

17. A method in accordance with Claim 1^ wherein 
the step of backprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes 
through }^'^ , 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
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tion data strip of the group of correction projec- 
tion data strips E? through E^^ in a pixel position 
influenced by the selected voxel multiplied by a 
weighting factor h(p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

18. A method in accordance with Claim 15, wherein 
the step of backpro jecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes Vo 
through V^'^ , 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data strip of the group of correction projec- 

<iata strips E? through E^'^ in a pixel position 
influenced by the selected voxel multiplied by the 
weighting factor h(p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

19. Parallel processing apparatus implementing the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 

5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
source scanning trajectory to obtain, at each of a plurality 
of discrete source positions 6/ on the source scanning tra- 
jectory, a 2D measured cone beam pixel data set P., said appa- 
10 ratus comprising: 

a data memory for storing the voxel data set V 
organized as a plurality m of independent voxel subcubes 
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through V^"^, each voxel subcube including a plurality of vox- 
els; 

15 processing elements for initializing the voxel data 

set V and performing successive iterations to correct the 
values of at least selected voxels of the 3D voxel data set 
based on the measured cone beam data sets P^, during each 

iteration correcting the values of at least the selected vox- 
20 els for each of the discrete source positions 0/ on the source 

scanning trajectory, said processing elements including in a 
level 0 a plurality m of backpro ject/repro jectors 
corresponding to the voxel subcubes, and in a level 1 at 
least one split/merge processor, 
25 said backpro ject/repro ject processors operable 

to reproject each voxel subcube through V^'^ to 
calculate pixel values for a corresponding 2D cal- 
culated projection data subset from a group of cal- 
culated projection data subsets P^^ through i^'""^ of a 
30 2D calculated projection data set Pi, the projec- 

tion data subsets P^ through P^'^ overlapping in 
part/ 

said at least one split/merge processor opera- 
ble to merge the calculated projection data subsets 
35 P^ through P{^'^ to calculate pixel values of the 2D 

calculated projection data set Pi, adding pixel 
values in any overlapping regions of the calculated 
projection data subsets P^ through /^'""^ 

said processing elements including means for 
40 calculating a 2D correction projection data set Ei 

by determining the normalized difference between 
the value of each pixel of the measured projection 
data set P^ and the value of the corresponding pixel 
of the calculated projection data set Pi, 

45 said split/merge processor operable to split 

the 2D correction projection data set £/ into a 
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plurality of 2D correction projection data subsets 
E- through E/""^ corresponding to the voxel subcubes 
vo through V^^-l, the 2D correction projection data 
50 subsets £; through E^' overlapping in part, dupli- 

cating element values for any overlapping regions 
of the correction projection data subsets E^ 
through E^^^ , and 

said backproject/reproject processors operable 
55 to backproject each correction projection data sub- 

set £f through ^ to correct voxel values in the 

corresponding voxel subcube of the plurality of 
voxel subcubes through V^'^ . 

20. Apparatus in accordance with Claim 19, which 
comprises a plurality of split/merge processors organized as 
a tree-type structure in a plurality of layers; 

said split/merge processors operable to merge 
5 groups of calculated projection data subsets from the next 
lower layer to present a lesser number of calculated projec- 
tion data subsets to the next higher layer, the lowest layer 
merging groups of the calculated projection data subsets 
through P^^"^ from the backproject/reproject processors, and 

10 the highest layer producing the calculated projection data 
set Pi; and 

said split/merge processors operable to split one 
or more correction projection data subsets from the next 
higher layer to present a greater number of correction pro- 
15 jection data subsets to the next lower layer, the highest 
layer splitting the correction data set £/, and the lowest 
layer producing the correction projection data subsets £? 
through £ . 

21. Apparatus in accordance with Claim 20, wherein 
at least said split/merge processors comprise transputers. 
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22. Apparatus in accordance with Claim 20, wherein 
said split/merge processors of each layer each include one 
I/O port connected to the next higher layer and three I/O 
ports connected to the next lower layer. 
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